#Plot of core and edge proportions (plus sds) with no parasite cost at k=32000 and cell.size =5

setwd("/Users/ben/Documents/Papers/Published/Rhabdias/Parasite lag model/Pathol Freq Ind Dyn Load")


y.error.bars<-function(x, y, ybar, lty="solid", ...){
	points(x, y, ...)
	arrows(x, y-ybar, x, y+ybar, code=3, angle=90, length=0.1, lty=lty, ...)
}

data<-read.table(file="Pathol Freq indepdt r15 Final.txt", header=T)
attach(data)

par(mfrow=c(2, 1)) #530 x 775
## First plot
core<-tapply(core.prop, list(surv.cost, incidence), mean)
edge<-tapply(edge.prop, list(surv.cost, incidence), mean)
core.sd<-tapply(core.prop, list(surv.cost, incidence), sd)
edge.sd<-tapply(edge.prop, list(surv.cost, incidence), sd)
x<-as.numeric(colnames(core))

par(mar=c(3,5,5,2))
plot(0:1, 0:1, ,type="n", bty="l", xlab="", mar=c(4,5,4,2), cex.lab=1.5, cex.main=1.5,ylab="Observed prevalence", main="Density independent transmisison")
cols<-c("black", "grey50", "grey30")
for (i in 1:3){
	if (i==1) y.error.bars(x, core[i,], core.sd[i,], lty="solid", pch=21, cex=2, col=cols[i])
	else points(x, core[i,], pch=19, col=cols[i], cex=2)
	lines(x, core[i,], col=cols[i])
}
for (i in 1:3){
	if (i==1) y.error.bars(x, edge[i,], edge.sd[i,], lty="solid", pch=24,cex=2, col=cols[i])
	else points(x, edge[i,], pch=17, col=cols[i], cex=2)
	lines(x, edge[i,], col=cols[i])
}

legend(0, 0.7, legend=c("Range core", "Range edge"), pch=c(21, 24),cex=1.3, bty="n")
legend(0, 1, legend=c("= 0", "= 0.1", "= 0.4"),fill=c(NA, cols[-1]), cex=1.3, bty="n", title="Cost of parasitism")

## Second plot

summary<-tapply(x.diff, list(surv.cost, incidence), mean)
summary.sd<-tapply(x.diff, list(surv.cost, incidence), sd)
x<-as.numeric(colnames(summary))
summary[which(summary==Inf)]<-NA
summary.sd[which(is.nan(summary.sd))]<-NA

lins<-c("solid", "dashed", "dotted")
par(mar=c(5,5,4,2))
plot(0:1, c(-4, max(summary, na.rm=T)), ,type="n", bty="l", cex.lab=1.5, xlab="Transmission rate", ylab="Lag distance (arbitrary units)")
for(i in 1:3){
	if (i==1) y.error.bars(x, summary[i,], summary.sd[i,], pch=21, cex=2)
	else points(x, summary[i,], pch=21, cex=2)
	lines(x, summary[i,], lty=lins[i])
}
legend(0.7, 90, legend=c("= 0", "= 0.1", "= 0.4"), lty=lins, cex=1.3, bty="n", title="Cost of parasitism")
## Clean up

detach(data)
rm(data, core, core.sd, edge, edge.sd, x, y.error.bars, summary, summary.sd)

